ez = 1;
ev = linspace(0,4,201)';

theta = linspace(0, 2*pi, 1001);

edepair = ez + ev*cos(2*theta);

eavg = mean(abs(edepair),2);
figure(9868);
plot(ev, eavg, ev, ev*2/pi, ':', ev, ez*ones(size(ev)), ':');
xlabel '\epsilon_V/\epsilon_Z'
ylabel '\epsilon_{depair}/\epsilon_Z'